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the basis of the extended Kalman filter (EKF) and the 
complimentary Kalman filter developed in Section 4.2. A 
discussion of Kalman filtering can be found in [12]. 

3. Literature analysis 

In applying Kalman filtering to the inertial orientation 
tracking problem there is considerable freedom in system 
modeling - what physical variables to assign to the state 
vector x, what measurements are in the measurement 
vector y, and what matrices A, B, C, Q, and R most 
accurately describe the system given those choices. A 
literature search was conducted to see how other authors 
have used Kalman filters to estimate orientation from the 
outputs of 3 strapdown gyros. The 7 most relevant 
references found are reviewed in this section. Two come 
from vehicle navigation, two from robotics, and three 
from virtual environments. 

An early maritime navigation work by Bona and Smay 

[13] , summarized in [12], is of interest because it shows 
how to reset gyro biases based on indirect measurements 
(position errors that result from them) and provides a 
now-common Markov model of gyro bias evolution. The 
dynamic system model details how the position errors 
evolve in response to the gyro biases, and how the gyro 
bias Markov components evolve in response to the 
process noise. 

The most relevant reference found in the aeronautics 
literature was Koifman and Merhav’s description of an 
autonomously aided strapdown attitude reference system 

[14] . Here, an autopilot is created with three low-cost rate 
gyros with time- varying biases on the order of 0.1 °/s. The 
measurements fed into the Kalman filter are from the 
three gyros, a magnetic compass, altimeter, and airspeed 
sensor. The state vector contains 16 elements: 3 linear 
velocities, 3 angular velocities, 3 orientation Euler angles, 
altitude, 3 wind gust velocity components and 3 gyro 
biases. The state transition matrix is obtained by lineariz- 
ing the system differential equations which encompass 
the aircraft equations of motion as well as the kinematic 
Euler equations (6). In contrast to Bona and Smay, the 
gyro biases are considered piecewise constant, and the 
corresponding diagonal covariance elements are simply 
reset whenever a change detection algorithm suspects that 
the gyro biases may have changed. It is also instructive to 
note that the full order 16-dimensional system could not 
be run in real time, so they reduced the state to 11 
elements and were then able to achieve about 20 updates 
per second with minimal loss in accuracy. The measure- 
ment vector consists of the three gyros and the airspeed 
sensor. 

Barshan and Durrant- Whyte [15] investigated the use 
of a solid-state gyroscope for mobile robotics applications. 


They paid particular attention to the gyroscope error 
model, and came up with an exponential curve to fit the 
changes in bias as the gyroscope warms up. They then 
implemented a Kalman filter for estimating a single 
rotation angle O, with a state vector 
d> £* and a state transition matrix 

that propagates the truth states and error states 

e <s> ,£ < p completely independently. The only system 

observation is the single rate gyro measurement, so the 
system is not observable, and the angular position error 
covariance grows unbounded. However, it is demonstrated 
that the gyro drift error grows at a rate 5 times slower 
when using the exponential gyro error model. 

A paper on mobile robot attitude estimation by Vaga- 
nay et al [16] provides the only example in the literature 
in which gyroscope drift is compensated using two 
accelerometers, and is therefore particularly germane to 
this drift-free head-tracking application. The Kalman 
filter model is very unusual and results in a state vector of 
surprisingly low dimension. The integration of angular 
rates is done outside of the Kalman filter, and is treated as 
part of a measurement system that provides gyroscopically 
determined measurements of pitch and roll, 0 g and y g , 
which are complimented by gravimetric measurements erf 
0 and \|f from the accelerometers. The state contains 0 and 
y and the pitch and roll drift rates, and the transition 
matrix used in the Kalman filter is just the identity. This 
is the leanest Kalman filter conceivable, as even the 
kinematics of Euler angle integration are not modeled, but 
the performance reported is nearly comparable to the 
other methods. No details are given about the determina- 
tion of Q and R. 

Azuma and Bishop developed a Kalman filter to use 
inertial sensors together with an optical head-tracker to 
predict head motion in HMD applications [17]. The 
approach is different from the preceding papers, and also 
from the application developed in this paper, because the 
gyroscope rate signals are not integrated to obtain 
orientation. Instead, the orientation is obtained from the 
optical head-tracker, and the angular rates are fused with 
this in the Kalman filter to yield improved predictions. 
The state vector contains a quaternion specifying 
orientation, the angular rates in body axes, and the 
angular accelerations in body axes. The measurement 
consists of the quaternion measured by the optical tracker, 
and the angular rates measured by the gyros. The Q and 
R matrices are determined off-line using Powell’s method 
on prerecorded datasets to find the parameters that give 
the best performance. Prediction was accomplished by 
extrapolating forward in time, using the angular velocity 
and acceleration estimates in the state vector. 
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Emura and Tachi likewise used gyros to augment the 
dynamic performance of an existing head-tracker, but in 
this case the tracker was magnetic instead of optical [18, 
19]. The state vector contains orientation (Euler angles in 
the first paper were replaced with a quaternion in the 
second) and angular velocities. The measurement vector 
measures all elements of the state, using a Polhemus 
magnetic tracker to measure orientation and gyros to 
measure the angular rates. A novel aspect of the Kalman 
filter structure is the use of two different types of meas- 
urement update step: a 3-dimensional measurement used 
most of the time, when only gyro data is available, and a 
6-dimensional measurement used when the Polhemus 
data is available as well. Q and R were found empirically, 
using a high-precision mechanical tracker as a reference 
to measure remnant errors. 

4. System modeling and filter design 

4.1 State and measurement vectors 

The first sten in modeline_is_to decide what tn_nut Jjv 


allows the gyro measurements to be utilized in the 
obvious way - as measurements. However, while it is 
obvious from (6) how the derivatives of the orientation 
state elements will be computed from the state, how shall 
the derivatives of the angular velocity components depend 
on state? Some authors [18, 19] simply assume zero 
dependence, i.e. constant angular rates. Some process 
noise is added to the angular accelerations to allow for 
non-constant angular rates, but in reality the angular 
accelerations would not be very much like white noises, so 
this model cannot be very optimal. Other authors [15, 17] 
augment the state vector with to, which changes the 
model to an assumption of constant angular acceleration. 
The difference between the true © and the assumed 
© = 0 is closer to white noise. Further derivatives, as in 
[15], make the model even more accurate, but lead to an 
unreasonably large state vector. 

For most accurate estimation, the equations of motion 
of the body being tracked should be included in the 
system dynamics model (1). For example, in [14] the 
angular accelerations of the aircraft depend precisely, 
through well-known aircraft equations of motion, on 
quantities in the state vector and aileron positions, which 
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Figure 3: Complimentary Kalman filter for orientation 

tees that the rapid dynamic response of the inertial system 
will not be compromised by the Kalman filter. Another 
advantage is that the gyro rates are not treated as 
measurements, so it is unnecessary to include 0) in the 
state vector. Since the head dynamics are not being 
modeled in this implementation, co is excess baggage, and 
by removing it from x the dimension is reduced from 9 to 
6, with more than a three-fold computational savings. The 
following sections, therefore, will strive to develop a 
complimentary Kalman filter to estimate 

^ X = [sco] S [^ ^ ^ ^ C ° x $ C ° y (7) 

using 

8y = inclinometer ~ W ® inclinometer ~ ® compass 

as the measurements, where 50 represents the error in the 
output of the attitude computer, and Sco represents the 
gyro biases. 

4.3 DT nonlinear attitude computation 

A Kalman filter which operates on the errors of the 
INS attitude computer must mimic the noise-free error 
dynamics of the attitude computation. This section derives 
the attitude integration algorithm, Section 4.4 linearizes 
the attitude algorithm to obtain the error dynamics, and 


Section 4.5 describes a complimentary 
EKF to operate on the errors of the attitude 
computation with the computational 
complexity of the EKF reduced by 

applying Friedland’s separate bias 

formulation. 

The continuous-time (CT) nonlinear 
& differential equation which the attitude 

^ computer must integrate was given in (6). 
To derive the DT attitude computation 
from it, it is useful to approximate the 
evolution of 0(t) over a short time interval 
by its Taylor series expansion 

A / 2 

0(/ + A/) = 0(r) + 0(r)A/ + 0(O— +... (9) 

The number of terms which must 
attitude 0 be retained depends on the size of 

m totinn I ^ At. For a first order integration 
"Om Py aTiQn J algorithm (retaining only the first 

§q two terms), the error per step will 

be mostly due to the third term, 

_ which is of order co 2 At 2 /2. 

timator <_ _ , 

Therefore, 

1 2 , 

error rate * —co A/. 

2 

For typical peak head velocities of 
about 6 radians/sec and a timestep 
of 0.003 sec, this yields an error rate of about .05 rad/s 
(about 3°/s) which is unacceptable. Retaining the third 
term, the error rate will be dominated by the fourth term, 
or order a) 3 Af 3 /6, so 

1 3a 2 

error rate «— co At . 

6 

For the same to and At the error rate would be about 
0.0003 rad/s, or about I°/min. Since the low-cost gyros 
are unlikely to have performance much better than this, a 
second order integration algorithm was selected. 

Differentiating (6) by the chain rule for partial deriva- 
tives results in 

®W“ 4 [w.( 0 (O) ®W] 9 (')+ 

(10) 

-£[w a (e(0) ©(*)] ©(0 

Defining (with time indices suppressed for brevity) 

v,(e.®)=|(w,(eH= 

cosyand sinyanO sny cosy 

Ct) a (0 + G) ( 0 

ccs0 * oos0 ccs 6 ax 6 (11) 

-any OJ -easy cd 0 0 

ccs y any sin y/ sin 0 cosy/ sin 0 

(o - co CD + cd 0 

ccs 6 ' cos 6 cos 6 cos $ 


sin y/ sin 0 cosy/ sin 0 
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and approximating the derivative of co(t) by its first 
difference. 


<0 (0- 


©(f + A/)-©(r) 
A t 


and substituting (1 1) and (12) into (10) yields 

0(O = v i (0(O.fi)(O) w a ( 6 (/)) ft )(0 


( 12 ) 


+w a (0(O) 


©(f + A/)-©(r) 
A/ 


(13) 


Plugging (6) and (13) into (9) and rearranging terms 
slightly leads to 


8(r + Ar) = 9(t) + W B — - ) ' + — + A/) Ar 


+V fl W B o)(r) 


A/ J 


(14) 


which is the second order DT integration step formula 
implemented in the attitude computer. Since At remains 
as an explicit parameter in this formula, it is unnecessary 
to have constant stepsize. This eliminates the difficulties 
of an interrupt driven program structure that would be 
necessary to have constant sampling rate data acquisition. 


4.4 DT linearized error dynamics 


Equation (14) defines a nonlinear state propagation 
function f At for the system with state vector 0 and input ©: 
8(f + A/) = f* (8(r), ©(f), ©(f + A/), f) (15) 


For the sake of obtaining an extended Kalman filter 
which can estimate both orientation errors and gyro 
biases, consider augmenting the state vector with © and 
rewriting the system in the form 


8(f 

©(f 


L ( 


+Ar)i ~ rr 6 (oT 

+ Af)J ^[©(f)J 


+ u(f) 


0(0 

©(f) 


(8(f), ©(/),©(/ + Af), Af) 

©(f) 


(16) 


“(0 = 


0 

to(r + At) - ©(/) 


where u(t) has been deviously chosen to make co(t) evolve 
in accordance with the input history of the previous 
system. The system error dynamics can now be obtained 
by linearizing about the nominal trajectory [0(r) a>(r)] r 


to get 


S8(/ + A/) TA 

. BT 88(/) 

8©(r + At) 8 

l]s©(/) 


where 


(17) 


A = 

B = 


jUO 

<*K0 

jUO 

<to(f) 


= W B Af+[v 8 W, + (^V a jw,© 


Af 2 

2 

A/ 2 

2 


(18) 


and 0 and I are 3-by-3 zero and identity matrices. The 
vector partial derivatives of V B are too messy to write out 
in full, but the computation is straightforward and can be 
carried out as follows: 1) form a “row vector** of the three 
matrices obtained by differentiating V B with respect to the 
first, second and third elements of the vector in the 
denominator of the partial derivative; 2) multiply each of 
these three matrices by the r.h.s. vector W B co. This results 
in a “row vector of column vectors**, i.e. a 3-by-3 matrix. 

Equation (17) gives the state transition matrix for the 
linearized error dynamics of the augmented system. The 
angular velocity errors 5© are principally due to gyro 
biases, and will be interpreted simply as gyro biases from 
here on. The A and B submatrices can be interpreted as 
describing the influence of the orientation error and gyro 
biases at time / on the orientation error at time t+At. The 
effect of the matrix is fairly obvious; it basically mimics 
the attitude computation of (14) except that the input 
angular velocity is due to gyro biases and the output is 
therefore an orientation error . The growth of orientation 
error in the absence of angular rate errors is governed by 
the A matrix. To first order A = I+V a A/. The identity 
term maintains the previously accrued error, and V B (0,<o) 
amplifies existing orientation errors in response to 
motion. 


4.5 Separate-bias Kalman filter formulation 

The linear error propagation model of (17) provides 
the basis for a complimentary Kalman filter to estimate 
these errors. The model has been manipulated into a form 
in which the gyro biases are assumed constant, thus 
permitting the direct application of Friedland’s separate- 
bias Kalman filtering results [21]. If the constant-bias 
model turns out to fit the gyro performance poorly, the 
restriction can later be ameliorated by use of an age- 
weighting factor. If an exponential gyro warm-up model 
as in [15] seems more appropriate, this can be accommo- 
dated within Friedland*s formulation by replacing the 
identity submatrix in the state transition matrix of (17). 

Switching to Friedland’s notation, define an error state 
vector x k h SQ(t k ) and a bias state vector b* s 5co(r A ) 
where t* is the time at the k* iteration of the algorithm. 
An augmented state vector z k = [x k bj r satisfies 
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The additive white noise w k , with variance Q k , only 
effects x, since b is assumed constant. The measurement 
equation is 

y t =L A z t +v t , (20) 

where v k is white noise with variance Q k . In Friedland’s 
paper, L t =[H t C t ], but in this application the 
measurements from the inclinometers and compass only 
measure x and not b, so C = 0 will be used throughout, 
resulting in a great simplification from Friedland’s 
derivation. 

Applying Kalman filtering to this model, the optimal 
estimate of z is 

£ t+I =F,z t + K(Jk+lXy t+l -LF t z,) (21) 

K (*) = P(fc)L r [LP(fc)L r + R* ]"' . (22) 

The Ricatti equations for the recursive computation of the 
estimation error covariance matrix P(k) needed in the 
Kalman gain expression can be rolled together into the 
single predictor-to-predictor covariance update equation: 

P(*+l) = F t [l-K(*)L]P(*)F/ +[J]Qw[I °]( 23 ) 

Partitioning P(k) into 3-by-3 submatrices as 

r p,<*> p.wi (24) 

() [p/W P.wJ’ 

the expression for the Kalman gain, (22), may be 
rewritten in partitioned form as 

[ P x (*)H r [HP,(*)H r + R 4 ]'‘ 1 

|p x /(*)H r [HP x (*)H r + R 4 ]~ J ‘ 

These separate gains are used in two essentially separate 
Kalman filters, one for estimating x and one for b. To 


2-axis fluid 
inclinometer 
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magnetometer 


compute the K x and K*> gains in (25), covariance sub- 
matrices P x and P xb are needed. These are updated by the 
partitioned version of (23): 

& Sl-R; ‘t :H3" •')■ 

rp x P^TA. B*] JQ. 0‘ 

P/ P, 1 0 I I I 0 0 

L - J L J (26) 

[A 4 -A 4 K x H-B 4 K t H B t ] 

"I "K fc H I J X 

[PA r +PX P ^1J Q ‘ °‘ 

|p j /A 1 r +P fc B 4 T P*J [0 0 

Thus, a plethora of 6-by-6 matrix multiplications and 
one 6-by-6 inversion are replaced by a somewhat greater 
number of 3-by-3 multiplications and one 3-by-3 inver- 
sion. 

5. Implementation 

Figure 4 illustrates the configuration of the hardware 
built to demonstrate the inertial head-attitude tracking 
concept The sensors are all embedded in a specially 
machined 2" X 2" X 1.25" plastic block connected by a 
thin 10' cable to an analog signal conditioning circuit and 
data acquisition card in a PC. 

Software was written in “C” to run on the PC and 
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Figure 4: Orientation tracker hardware configuration. 


Figure 5: Inertial orientation tracker main software loop. 

The initialization block, executed once 
at program start-up, sets the initial state 
estimates and covariances as follows: 

Xo: The inclinometer is read and used 

6dx33 to set y q com pa 5 Si if used, 

i r T 

determines <{>; otherwise 4>=0- 

ie |, bo- The biases of all 3 gyros are meas- 

\/D ured during system calibration and stored 

ter card in a file. On initialization, the file is read 

and 5oo is initialized with the stored gyro 

biases. 

P x (0): The errors in the initial deter- 
mination of the Euler angles may be 
substantial, but they are assumed to be 
uncorrelated with one another: P x (0) = I. 
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P b (0): The gyro biases at start-up could differ substan- 
tially from the prerecorded calibration values, but the 
uncertainties are uncorrelated: P b (0) = 0.1 1. 

P*b(0): The initial uncertainties in orientation and gyro 
bias are completely uncorrelated: P xb (0) = 0. 

The data acquisition block scans all the A/D channels 
in rapid succession. The new gyro readings are stored as 
©(t+At) and the previous ones are moved back to ©(t). 
The new inclinometer and compass readings are stored in 
y(t+At). In the next block, a timestamp is obtained from 
the 8253 timer/counter chip on the PC motherboard. This 
counter is driven by a 1.19 MHz oscillator with a 65,536 
divisor to generate 18.2 Hz timer ticks for BIOS and DOS 
time-keeping. By reprogramming the divisor it was found 
possible to obtain sub- microsecond timing resolution as 
required for inertial integration. At is calculated as the 
difference between the current timestamp and the 
previous one. 

Next, co(t), ©(t+At) and At are fed into the Kalman 
filter update block. W B and V B are computed and then 
used in (14) to compute the predicted 0(t+At). This 
corresponds to the attitude computation block in. Since 
the Euler angle estimates, 0 must be maintained anyway, 
it is convenient to subsume 50 into them, and keep track 
of total estimates only. This does not change the filter 
framework developed in the previous section in any 
important way; it just means that 80(/) is always zero at 
the beginning of each iteration of the Kalman filter. At 
the end of the Kalman filter update cycle, 50(f + Af) is 
used to reset 0(f + Ar) and then flushed back to zero 
before the next cycle. Since the attitude error estimates 
are propagated along with the attitude estimates through 
the nonlinear propagation equation, the top three 
elements of F k z k in (21) are replaced with zeros. Since © 
is not included in the state, the running estimates of 8© 


covariance submatrices using (26). Since the inclinometer 
and compass signals are pre-processed to give direct 
measurements of the Euler angles, H=I, and (26) is 
simplified to the following steps: 

T, = A-AK, 

t 2 =t,p* 

T 3 = BT/ 

P,* = Tj + BP, (28) 

P^=T J+ BP/ 

P* = P lt + B r +T 3 +T,P x A r 

where Tj are simply temporary storage matrices used to 
reduce the amount of redundant matrix multiplication. A 
small subroutine library was written, following the pointer 
conventions and numerical methods described in [22], to 
perform the necessary matrix multiplication, transposi- 
tion, addition and inversion operations to carry out these 
steps. 

5.1 The Q k and R k Matrices 

Ideally, Q k is supposed to reflect the magnitude of a 
white noise sequence. If all error sources in the inertial 
attitude system are taken care of (i.e. modeled in the state 
propagation matrix), then w k in (19) should be entirely 
due to the noise floors of the angular rate sensors. In this 
case, it should be possible to calculate the optimal value of 
Q k by measuring the noise covariance, Q, of the station- 
ary gyros in advance, then at each time step compute 
Q* = G.QG/, using G* = W,(0(f t )). 

However, there are many nonwhite error sources be- 
sides bias, such as nonlinearity, hysteresis, misalignment, 
g-sensitivity, and scale factor temperature coefficient, 
none of which are modeled in the current implementation. 
The best procedure for designing a reduced-order Kalman 
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scale factor error: This is a composite of nonlinearity 
and temperature dependent scale factor variations. 
Assuming scale factor accuracy of 1% of full scale, 
a=0.01co rad/s. Covariance per step (O.OlcoAt) 2 . 

Assuming At=0.01sec, and that these error sources are 
uncorrelated, the error covariances add up to approxi- 
mately lO'^l+oVlO* 4 ® 6 ). At each iteration of the 
Kalman filter software, the following algorithm is used to 


4. 


set R k 





0 0 

°r 7 0 . 

0 a w 


According to this algorithm, the measurement error 
covariances for the inclinometer roll and pitch range from 
1, during periods of likely slosh, down to 10 -4 , during 
periods of likely stillness. The covariance of the compass 
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Kalman filter block is disabled by setting K x and Kb equal 
to zero. During the test period of approximately 35 
seconds, the sensor block was repeatedly turned through 
+90° about the roll axis and left to rest on its right side, 
then returned to rest in its horizontal orientation on the 
table. The roll Euler angle is plotted against time in 
Figure 6, which demonstrates the problem with unaided 
inertial integration: the accumulated drift error by the end 
of the run is about 15°. The second dataset is created by a 
similar motion sequence, but the Kalman filter is in 
effect As Figure 7 shows, the filter incorporates the drift- 
free but noisy measurements from the inclinometers, and 


[4] K. R. Britting, Inertial Navigation Systems Analysis. 
New York: Wiley-Interscience, 1971. 

[5J C. Broxmeyer, Inertial Navigation Systems. New 
York: McGraw-Hill, 1964. 

[6] J. L. Farrell, Integrated Aircraft Navigation. New 
York: Academic Press, 1976. 

[7] A. Lawrence, Modem Inertial Technology: Springer- 
Verlag, 1993. 

[8] R. H. Parvin, Inertial Navigation. Princeton, New 
Jersey: Van Nostrand, 1962. 

[9] G. M. Siouris, Aerospace Avionics Systems. A Mod- 
em Synthesis. San Diego, CA: Academic Press, 1993. 

[10] E. Foxlin, "Inertial Head-Tracking," M.S. Thesis, 
Dept of Elec. Engineering^ and Cbyng^ yj^ M^ss # 


Proceedings of the 

* 


IEEE 1996 Virtual Reality Annual 

International Symposium 


March 30 - April 3, 1 996 Santa Clara, California 


Sponsored by 

The IEEE Computer Society Technical Committee on Computer Graphics 
The IEEE Neural Networks Council Virtual Reality Technical Committe 


IEEE Computer Society Press 
Los Alamitos, California 


Washington 


Brussels 


Tokyo 




IEEE Computer Society Press 
10662 Los Vaqueros Circle 
P.O. Box 3014 

Los Alamitos, C A 90720-1264 


Copyright © 1996 by The Institute of Electrical and Electronics Engineers, Inc. 

All rights reserved. 


Copyright and Reprint Permissions : Abstracting is permitted with credit to the source. Libraries may 
photocopy beyond the limits of US copyright law, for private use of patrons, those articles in this volume that 
carry a code at the bottom of the first page, provided that the per-copy fee indicated in the code is paid 
through the Copyright Clearance Center, 222 Rosewood Drive, Danvers, MA 01923. 

Other copying, reprint, or republication requests should be addressed to: IEEE Copyrights Manager, IEEE 
Service Center, 445 Hoes Lane, P.O. Box 1331, Piscataway, NJ 08855- 1331. 

The papers in this book comprise the proceedings of the meeting mentioned on the cover and title page. They 
reflect the authors’ opinions and, in the interests of timely dissemination, are published as presented and 
without change. Their inclusion in this publication does not necessarily constitute endorsement by the 
editors , the IEEE Computer Society Press, or the Institute of Electrical and Electronics Engineers. Inc. 


IEEE Computer Society Press Order Number PR07295 
ISBN 0-8186-7295-1 
ISSN 1087-8270 

IEEE Order Plan Catalog Number 96CB35922 
IEEE Order Plan ISBN 0-8186-7296-X 
Microfiche ISBN 0-8 1 86-7297-8 


Additional copies may be ordered from: 


IEEE Computer Society Press 
Customer Service Center 
10662 Los Vaqueros Circle 
P O. Box 3014 

Los Alamitos, CA 90720-1264 
Tel: +1-714*821-8380 
Fax: +1-714-821-4641 
Email: cs.books@computer.or 


IEEE Service Center 
445 Hoes Lane 
P.O Box 1331 
Piscataway, NJ 08855-1331 
Tel: +1 -908-981-1393 
Fax: +1-908-981-9667 
mis.custserv@compuier.org 


IEEE Computer Society 

13, Avenue de PAquilon 

B- 1 200 Brussels 

BELGIUM 

Tel: +32-2-770-2198 

Fax: +32-2-770-8505 

curo.ofc@cpmputer.org 


IEEE Computer Society 
Ooshima Building 
2-19-1 Minami-Aoyama 
Minato-ku, Tokyo 107 
JAPAN 

Tel: +81-3-3408-3118 
Fax: +8 1-3 3408-3553 
tokyo.ofc ©computer org 


Editorial and cover production by Mary E. (Cavanaugh 
Printed in the United States of America by KNI Inc. 


Y 

:,)Ci 



